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We  test  the  behavior  of  a  unified  continuous/discontinuous  Galerkin 
(CG/DG)  shallow  water  model  in  spherical  geometry  with  curved  elements 
on  three  different  grids  of  ubiquitous  use  in  atmospheric  modeling: 
(A)  the  cubed-sphere,  (B)  the  reduced  latitude-longitude,  and  (C)  the 
icosahedral  grid.  Both  conforming  and  non-conforming  grids  are  adopted 
including  static  and  dynamically  adaptive  grids  for  a  total  of  twelve  mesh 
configurations.  The  behavior  of  CG  and  DG  on  the  different  grids  are 
compared  for  a  non-linear  mid-latitude  perturbed  jet  and  for  a  linear  case 
that  admits  an  analytic  solution.  Because  the  inviscid  solution  on  certain 
grids  shows  a  high  sensitivity  to  the  resolution,  the  viscous  counterpart 
of  the  governing  equations  are  also  solved  and  the  results  compared.  The 
logically  unstructured  element-based  CG/DG  model  described  in  this  paper 
is  flexible  with  respect  to  arbitrary  grids.  However,  we  were  unable  to  define 
a  best  grid  configuration  that  could  possibly  minimize  the  error  regardless 
of  the  characteristic  geometry  of  the  flow.  This  is  especially  true  if  the 
governing  equations  are  not  regularized  by  the  addition  of  a  sufficiently  large, 
fully  artificial,  diffusion  mechanism,  as  will  be  shown.  The  main  novelty  of 
this  study  lies  in  the  unified  implementation  of  two  element-based  Galerkin 
methods  that  share  the  same  numerical  machinery  and  that  do  not  rely  on 
any  specific  grid  configuration  to  solve  the  shallow  water  equation  on  the 
sphere. 
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1.  Introduction 

As  computational  power  increases,  meteorologists  demand 
greater  resolution  from  global  atmospheric  models  utilizing 
spherical  domains.  Driven  by  the  need  to  select  the 
most  proper  computational  grid  on  the  sphere  for  the 
Nonhydrostatic  Unified  Model  of  the  Atmosphere  (NUMA) 
(Kelly  and  Giraldo  2012;  Giraldo  et  al.  2013),  in  this  paper 
we  analyze  the  continuous  Galerkin  (CG  or  spectral 
element(s),  from  now  on)  and  discontinuous  Galerkin  (DG) 
solutions  of  the  shallow  water  equations  (SWE)  on  a  set  of 
common  grids  used  in  global  circulation  models  (GCMs). 
Specifically,  we  study  1)  the  cubed-sphere  (Sadourny  1972), 

2)  the  icosahedral  grid  (Ico)  (Sadourny  et  al.  1968),  and 

3)  a  type  of  reduced  latitude-longitude  (Lat-Lon)  geometry 
(Phillips  1957).  Because  SWE  capture  many  of  the  essential 
features  of  GCMs  while  eliminating  the  vertical  structure  of 
the  atmosphere,  the  results  from  this  study  will  be  directly 
applicable  to  the  fully  three-dimensional  model  NUMA. 

The  geometry  of  the  cubed-sphere  and  Lat-Lon  grids 
can  be  advantageous  for  a  finite  difference  based  solver, 
and  are  generally  usable  by  grid  point  methods  such 
as  finite  and  spectral  element  (FE,  SE),  DG,  or  finite 
volume  (FV)  methods.  Furthermore,  Lat-Lon  is  the  standard 
discretization  for  GCMs  based  on  spectral  transform 
methods  due  to  the  fast  Fourier  transform  operation  along 
the  longitude  direction  (Hogan  et  al.  1991).  Icosahedral 
grids  or,  more  generally,  unstructured  tri-  and  quad-based 
tessellations* *  of  the  sphere  are  geometrically  flexible,  but  not 
all  numerical  methods  are  able  to  handle  them.  Recently, 
different  grids  to  be  used  with  finite  volume  discretization 
techniques  have  been  analyzed  in,  e.g.,  Qaddouri  et  al. 
(2012)  or  Peixoto  and  Barros  (2013).  In  the  context  of 
Galerkin  methods,  however,  not  much  analysis  has  been 
conducted  with  respect  to  the  solution  dependence  on  these 
computational  grids.  Rather,  given  a  specific  grid,  we  are 
likely  to  find  a  model  that  is  developed  around  it  and  is 

'  Please  ensure  that  you  use  the  most  up  to  date  class  file,  available  from 
the  QJRMS  Home  Page  at 

http://onlinelibraxy.wiley.com/journal/10. 1002/ (ISSN) 1477-870X 

*The  keywords  grid,  mesh ,  and  tessellation  will  be  used  interchangeably 

throughout  the  paper. 
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optimized  to  minimize  the  error  in  a  specific  configuration. 
In  the  case  of  the  solution  of  shallow  water  problems, 
examples  of  this  approach  are  the  works  by  Taylor  et  al. 
(1997)  (SE  on  cubed-sphere),  Giraldo  (2001)  (SE  on  the 
quad-based  icosahedral  grid),  Giraldo  et  al.  (2002)  (DG  on 
a  triangle-based  icosahedral  grid),  or  Nair  et  al.  (2007)  (DG 
on  the  cubed  sphere).  Recently,  a  comparison  between 
dynamically  adaptive  and  uniform  meshes  was  performed 
by  Muller  et  al.  (2013)  for  triangle  based  discontinuous 
Galerkin  methods.  A  grid  comparison  using  spectral  elements 
and/or  finite  volumes  are  found  in  St-Cyr  et  al.  (2008)  and 
Weller  et  al.  (2009).  However,  the  current  work  presents 
relevant  differences  in  scope  and  goal  with  respect  to  these 
two  articles.  The  differences  can  be  listed  as  follows:  St- 
Cyr  et  al.  compare  two  different  numerical  methods  that, 
by  construction,  are  defined  on  two  different  grids.  Their 
comparative  analysis  was  made  in  terms  of  how  a  spectral 
element  shallow  water  code  and  a  finite  volume  code  compare 
with  each  other.  Their  main  concern  was  the  verification 
of  two  codes  when  adaptive  mesh  refinement  (AMR)  was 
used  with  no  emphasis,  in  that  work,  on  the  flexibility  of 
either  spectral  elements  or  finite  volumes  with  respect  to 
the  grid  structure.  On  the  other  hand,  Weller  et  al.  (2009) 
show  a  finite  volume  analogous  of  what  we  partially  show 
iu  this  work  using  spectral  elements.  Since  the  main  interest 
of  our  ongoing  work  in  global  atmospheric  modeling  is  that 
of  using  Galerkin  methods  that  should  not,  by  design,  be 
tied  to  any  grid  so  that  any  optimized  mesh  generation  tool 
could  be  used  to  build  our  grids,  we  found  that  we  needed 
a  comparison  along  the  same  lines  of  Weller  et  al.,  but 
taken  from  the  point  of  view  of  continuous  and  discontinuous 
Galerkin  methods. 

Our  study  is  performed  by  solving  the  two  non-linear  zonal 
flows  proposed  by  Galewsky  et  al.  (2004)  first,  followed  by 
the  test  case  2  by  Williamson  et  al.  (1992)  which  admits  the 
analytic  solution  of  a  steady  state  nonlinear  geostrophic  flow. 

In  the  case  of  the  Galewsky  flow  at  low  resolution,  and 
because  of  the  misalignment  of  the  flows  on  the  cubed- 
sphere  and  icosahedral  grid,  the  solution  diverges  from  the 
true  solution  in  the  case  of  (i)  the  equilibrium,  and  (ii)  the 
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barotropic  instability  tests.  As  the  resolution  is  increased,  the 
effect  of  the  irregular  grid  geometry  is  partially  camouflaged 
and  eventually  disappears.  In  search  of  the  coarsest  allowable 
resolution  to  be  used  while  still  preserving  accuracy,  we 
add  an  artificial  diffusion  with  constant  coefficient  v  = 
1  x  105nr2s-1  as  suggested  by  Galewsky  et  al.  (2004).  The 
viscous  simulations  are  defined  for  CG  only;  no  viscosity  is 
used  for  DG. 

Finally,  to  assess  the  ability  of  the  model  to  handle 
conforming  and  non-conforming  grids  with  refinement,  we 
complete  the  study  by  building  six  statically  adapted  grids. 
We  construct  a  fixed  high-resolution  grid  core  in  the  region 
where  the  zonal  jet  is  most  likely  to  develop,  and  coarsened 
the  rest  of  the  domain  with  a  1-  and  2-level  de-refinement 
approach  to  estimate  the  error  when  the  total  number 
of  grid  points  is  drastically  reduced.  The  DG  solution 
algorithm  on  the  refined  grids  is  based  on  the  procedure  by 
Kopera  and  Giraldo  (2014)  and  references  therein. 

The  paper  is  organized  as  follows.  In  Sec.  2,  we  report  on 
the  construction  of  high-order  grids  on  the  sphere.  The  model 
equations  are  described  in  Sec.  3  whereas  the  numerical 
method  of  solution  is  described  in  Sec.  4.  We  describe  the 
tests  and  analyze  the  results  in  Sec.  5.  Finally,  we  summarize 
our  findings  in  Sec.  6. 

2.  High  order  grid  generation  on  the  sphere 

Grid  generation  on  the  sphere  is  a  relatively  easy  task. 
However,  finding  the  grid  that  can  suit  different  numerical 
methods  is  not.  Because  element-based  Galerkin  methods 
have  the  advantage  of  being  flexible  with  respect  to  the  grid, 
we  analyze  how  different  meshes  on  the  sphere  can  affect  the 
spectral  element  and  discontinuous  Galerkin  solution  of  the 
shallow  water  equations. 

A  standard  approach  to  high  order  grid  generation  is  based 
on  the  construction  of  a  linear  grid  first  (i.e.  a  grid  with 
straight  edges  whose  extrema  are  the  only  points  that  lie 
on  the  geometry),  and  then  populate  it  with  the  high  order 
internal  points.  The  population  step  is  performed  element¬ 
wise  in  the  logical  space  obtained  by  a  proper  projection 
from  the  physical  space.  Once  the  higher  order  grid  points 
(c)  2013  Royal  Meteorological  Society 


have  been  built  on  the  plane,  a  backward  projection  onto  the 
physical  geometry  moves  the  new  high  order  elements  onto 
the  sphere.  This  process  is  such  that  the  high  order  elements 
approximate  the  surface  faithfully.  The  projection  from  the 
sphere  to  the  logical  space  may  differ  from  grid  to  grid.  The 
gnonronic  projection  by  Ronchi  et  al.  (1996)  is  used  to  build 
the  high  order  grids.  In  what  follows,  we  describe  how  the 
reduced  Lat-Lon,  cubed-sphere,  and  icosahedral  grids  are 
constructed. 

2.1.  Reduced  Lat-Lon  grid  (RLL) 

Since  the  effort  of  Phillips  (1957)  to  reduce  the  singularity 
problem  at  the  poles  of  a  global  latitude-longitude  grid, 
different  types  of  reduced  Lat-Lon  grids  have  been  used. 
Partially  based  on  the  composite  methods  of  Starius  (1977) 
and  Browning  et  al.  (1989),  Lanser  et  al.  (2000)  solved  the 
shallow  water  equations  on  a  Lat-Lon  grid  away  from  the 
poles  combined  with  a  stereographic  grid  at  the  poles.  In 
this  paper,  we  build  and  test  a  high  order  version  of  Lanser’s 
grids.  To  build  the  RLL  grid  we  use  a  multiblock  approach 
(see,  e.g.,  Thompson  (1987))  that  consists  of  building  a  set 
of  independent,  simply  connected  surface  grids  that  are 
eventually  patched  together  to  form  the  global  mesh.  The 
main  RLL  region  and  the  polar  caps  are  the  first  blocks 
to  be  built.  The  interfaces  between  the  RLL  area  and  the 
polar  patches  are  obtained  by  a  transfinite  interpolation 
(TFI,  Gordon  and  Hall  (1973);  Eriksson  (1984)).  The  way 
this  is  done  will  be  described  shortly.  The  main  Lat-Lon 
region  is  composed  of  four  faces  obtained  from  the  master 
face  71  =  [-71-/4  <  A  <  7t/4]  x  <  ip  <  pmax\  and  from 

its  rotation  with  respect  to  the  z-axis  of  the  sphere.  The 
variables  A  and  p  indicate  the  longitude  and  latitude, 
respectively.  The  RLL  band  in  the  northern  and  southern 
hemispheres  is  delimited  by  pmin  and  Pmax-  The  rotation 
matrix  from  71  to  the  three  remaining  faces  72,73,74  is 
obtained  by  the  combined  effect  of  a  translation  on  the  xy- 
plane  and  a  rotation  around  the  z-axis;  this  transformation 
is  given  by 
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The  first  matrix  produces  a  translation  [Imn]  =  [— 1  —  1  — 
1]  along  ( x ,  y,  z)  and  is  followed  by  the  rotations  \rot,i=2,z,4,  = 
(■7r/2,  7t,  37t/2)  given  by  the  action  of  the  second  matrix.  The 
coordinates  on  each  rotated  planar  face  are  simply  a  function 
of  ( x ,  y ,  2)7l  and  are  computed  as 

[x  y  Z  l]7i  =  [x  y  z  l]7l[T],  i  =  2,3, 4.  (1) 

Once  rotated,  each  gridded  face  is  projected  onto  the 
sphere.  The  polar  caps  are  built  in  the  same  way,  except 
that  the  starting  point  is  the  square  region  7 poiar  =  [— 7t/4  < 
A  <  7t/4]  x  [— 7t/4  <  ip  <  7r/4].  At  this  point  we  have  a  linear 
grid  made  of  6  disjoint  regular  patches.  The  top  view  of  this 
is  shown  in  Fig.  1.  Eight  new  surfaces  must  now  be  built  to 
fill  the  un-gridded  gaps. 


Figure  1.  Linear  grid.  Top  x  —  y  view  of  the  disjoint  lateral  Lat-Lon 
and  top  patches  in  the  multi-block  RLL  grid.  TFI  is  used  to  connect 
them  together  through  an  interface  surface  grid. 

To  build  each  interface  surface  by  TFI  we  need  information 
from  its  four  boundary  edges.  With  reference  to  the 
schematic  of  Fig.  2,  the  region  delimited  by  edges  1,2,  3, 4 
is  built  in  two  steps:  (a)  given  the  (A ,tp)  coordinates  of  the 
corner  points  1-2-3-4,  build  edges  3  and  4  and  subdivide  them 
into  Neiat  int  ID  linear  elements,  where  Neiat  int  is  the  user- 
defined  number  of  elements  along  the  latitude  direction  in  the 
©2013  Royal  Meteorological  Society 


interface  region;  (b)  compute  the  grid  points  in  the  interior  of 
the  patch  using  the  planar  and  linear  transfinite  interpolation 
defined  by  the  Boolean  sum 

X(|,t))=  U  +  V  — U®V,  (2) 

where,  given  the  arrays  £  =  [0, . . . ,  1]  and  rj  =  [0, . . . ,  1] 
in  computational  space  along  the  local  (u,  v)  directions, 
the  univariate  interpolations  and  tensor  products  (®)  are 
computed  as 

U(|,  fj)  =  (1  —  £)X(0,  v)  +  £X(1,  fj) 
V(|,f))  =  (l-r))X(?,0)+rjX(|,l) 

U  ®  V(|,  f})  =  (1  -  |)(1  -  77)X(0, 0)  +  (1  -  |)r)X( 0, 1) 

-(l-r))CX(l,0)-|l7X(l,l).  (3) 

In  Equations  (2),  X  is  the  array  of  the  (A,  95)  coordinates 
of  the  grid  points  in  the  interior  of  the  surface, 
given  the  known  values  of  the  surface  boundary  edges 
stored  in  X(0,  fj),  X(l,  77) ,  X(£,  0),  X(£,  1),  and  corners 
X(0, 0),  X(l,  1),  X(l,  0),  X(0, 1). 


The  high-order  Legendre-Gauss-Lobatto  (LGL)  points  are 
added  to  the  linear  grid  by  projecting  the  linear  elements 
onto  the  auxiliary  gnomonic  space  (on  the  plane),  and  back- 
projecting  the  new  high  order  quads  onto  the  sphere.  The 
total  number  of  elements  Ne  and  grid  points  Np  of  the  high 
order  conforming  grid  are  given  by: 
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AN eionN  i^N eiatN  -(- 1)  Lateral 

+2  {NeionN  +  l)(NeionN  +  1)  Polar 

+2(ANeionN(Neiat  intN  +  1)  —  8NeioriN)  Transition. 
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Figure  3.  Conforming  RLL  grid. 


In  (4)  and  (5),  Neiori  is  the  number  of  elements  of  order  N 
along  the  longitude  direction  of  one  face  in  the  main  Lat- 
Lon  band  that,  by  construction,  coincides  with  the  number 
of  longitude  elements  in  the  interface  and  polar  regions. 
Neiat  and  Neiai  int  are,  respectively,  the  number  of  latitude 
elements  of  the  Lat-Lon  region  and  of  the  interface  patches. 
The  number  of  latitude  elements  in  the  polar  caps  is  also 
given  by  Ne lon.  Fig.  3  shows  an  example  of  a  high-order 
conforming  RLL  tessellation. 

We  want  to  point  out  that  this  procedure  to  build  the 
multi-block  spherical  grid  is  not  unique  and  can  be  certainly 
optimized  to  give  the  best  distribution  of  the  high  order 
grid  points.  The  goal  of  our  study  is  not  to  propose  the 
most  optimized  grid  generators  but  rather  to  analyze  how  an 
element-based  Galerkin  model  handles  arbitrarily  generated 
grids. 

2.2.  Cubed- sphere 


LGL  points  are  omitted  for  the  sake  of  a  clearer  plot.  By 
construction,  the  elements  that  are  farthest  from  the  center 
of  each  face  are  smaller  than  those  at  the  face  center;  in 
addition,  their  distorted  shape  prove  to  be  a  concern  when 
the  resolution  is  not  sufficiently  high.  We  will  discuss  this 
issue  shortly. 


The  cubed  sphere  (see,  e.g.  Ronchi  et  al.  (1996);  Taylor  et  al. 
(1997,  2013))  is  derived  from  the  projection  of  a  cube  onto  the 
sphere.  As  such,  it  consists  of  6  faces  that  are  then  subdivided 
into  as  many  elements  as  necessary.  Like  RLL,  we  built  this 
grid  in  a  multi-block  fashion  by  going  through  the  same 
steps  described  above.  The  fundamental  difference  is  that 
the  cubed-sphere  grid  (also  referred  to  as  Hex  in  the  Figures 
and  Tables  throughout  the  paper)  has  only  6  faces  and  they 
are  all  equal.  An  example  is  shown  in  Fig.  4.  The  internal 
©2013  Royal  Meteorological  Society 


2.3.  Quad-based  icosahedral  grid 

The  quad-based  icosahedral  grid  of  order  N  is  derived  from 
an  initial  icosahedron  of  20  triangles.  For  better  properties  of 
the  high  order  grid,  every  triangle  is  mapped  onto  a  gnomonic 
space  by  the  transformation  formulas  (Giraldo  (2001)) 

*  =  ?’tanA/,  (6a) 
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y  =  r  tan  ip  sec  X' , 


(6b)  3.  Shallow  water  equations  on  the  sphere 


where  r  is  the  radius  of  the  sphere  and,  given  the  centroids  31  Equations  of  motion 
(Ac,  pc)  of  each  triangle,  we  define 


A'  =  atan 


cos</?sin(A  —  Ac) 


cos  ipc  sin  tp  +  cos  ipc  cos  ip  cos(A  —  Ac) 


(7a) 


p'  =  asin  [cos  ipc  sin  ip  —  sin  ipc  cos  p  cos(A  —  Ac)] .  (7b) 


After  mapping,  the  triangles  are  subdivided  into  smaller 
ones  by  a  Lagrange  polynomial  of  order  nj. 

The  number  of  quadrilateral  elements  and  grid  points  of 
the  final  grid  are  then  given  by 


Under  the  assumption  of  a  shallow  depth,  the  motion  of 
an  incompressible  fluid  may  be  described  by  the  non-linear 
shallow  water  equations.  The  shallow  water  assumption  holds 
as  long  as  the  characteristic  horizontal  extension  of  the  fluid 
is  much  larger  than  its  depth.  The  vector  representation 
of  the  shallow-water  equations  in  the  Cartesian  coordinates 
x  =  ( x ,  y ,  z)  is  given,  in  flux  form,  by 

^  +  V  •  (0u)  =  0  (8a) 

— f  V  ■  (cpu  ®  u)  =  — 0V0  —  /  (x  x  <£u)  +  px  +  5vS/2{4> u) 

(8b) 


Np  =  6 (IVj  -  2)N2  +  2 
Ne  =  6 (N%  -  2) 

where  the  number  of  points  NT  of  the  original  triangular 
grid  can  be  found  in  Giraldo  (2001). 

Fig.  5  shows  the  8th  order  icosahedral  grid.  As  for  the 
previous  cases,  the  elements  are  curved  and  lie  on  the 
spherical  surface. 


where  cf>  =  gh  is  the  geopotential  height  (g  and  h  are  the 
modulus  of  the  acceleration  of  gravity  and  the  vertical 
height  of  the  fluid),  vV2  is  the  artificial  viscosity  term  of 
viscous  coefficient  v  =  1  x  105  nr2  s-1  that  is  de-activated  by 
the  switch  5  =  0,  u  =  (u,v,w)'r  is  the  velocity  vector,  V  = 
( dx ,  dv,  dzy ,  and  /  =  2c oz/r2  is  the  Coriolis  parameter  with 
Earth’s  angular  velocity  uj  =  7.292  x  10-5  s~x  and  mean 
radius  r  =  6.37122  x  106m.  For  the  fluid  particles  to  remain 
on  the  spherical  surface  defined  in  3D  Cartesian  space, 
the  normal  component  of  the  velocity  field  with  respect  to 
the  sphere  is  removed  by  the  term  px  in  the  momentum 
equation,  where  p  is  the  Lagrange  multiplier  that  serves 
this  purpose.  The  construction  of  the  Lagrange  multiplier 
is  described  in  Sec.  3.2. 

3.2.  Forcing  the  fluid  onto  the  sphere  by  a  Lagrange 
multiplier 

We  use  the  approach  of  Cote  (1988).  Specifically,  at  the  end 
of  every  time  step  n  +  1,  the  momentum  of  the  flow  must  be 
corrected  for  the  particles  to  remain  on  the  sphere.  If  c  and  u 
denote  the  constrained  and  unconstrained  values  of  4>u,  the 
new  momentum  is  corrected  according  to 


(<H"+1  =  (<H"+i  +  a*- 


n+1 


(9) 
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For  the  fluid  to  remain  on  the  sphere,  the  condition 
u  •  x  =  0  must  hold  (i.e.,  the  velocity  field  and  the  position 
vector  must  be  orthogonal) ,  hence  implying  that  <puc  ■  x  =  0 
in  Equation  (9).  For  this  to  be  true,  the  value  of  the  multiplier 
is  found  to  be  p,  =  —<j> u"+1  •  x/r2.  Substituting  this  into  Ecp 
(9)  yields 


(<+) 


n+1 

c 


V  {<t> u 


n+1 


(10) 


where 


V  =  l- 


xxr 

r2 


(ii) 


is  an  orthogonal  projector.  T  indicates  the  transpose 
operator. 

The  Galerkin  discretization  of  (8)  is  described  in  the  next 
section. 


4.  Numerical  method 


4-1-  Spectral  element  and  discontinuous  Galerkin 
approximations 

To  solve  the  shallow  water  equations  by  element-based 
Galerkin  methods  on  a  domain  fi,  we  proceed  by  defining 
the  weak  form  of  (8)  that  we  first  write  in  compact  form  as 


^  +  V-F(q)  =  S(q),  (12) 

where  q  =  [<p,  <p u+  is  the  array  of  the  solution  variables,  and 
F  and  S  are  the  flux  and  source  terms  that  can  be  easily 
identified  in  the  equations  (8)  above.  Given  a  basis  function 
ip  that  belongs  to  the  Sobolev  space  H 1  in  the  case  of  spectral 
elements,  and  to  L 2  in  the  case  of  DG,  the  weak  form  of  (12) 
is  given  by  the  integral 


ip 


dq 

~dt 


+  V  ■  F(q)  —  S(q) 


dft  =  0, 


(13) 


that,  after  integration  by  parts  of  the  flux,  becomes 


f  ip-rr  dSl  —  j  ■  F(q)  dfl  —  f  ipS(q)  dfi  + 
Jn  Jn  Jn 


ipn  ■  F(q)  dT  =  0,  (14) 


where  F(q)  indicates  the  numerical  flux  of  the  advection 
and  diffusion  terms  across  the  element  boundary  T  of  unit 
outward  normal  n.  In  the  case  of  spectral  elements,  where  the 
basis  functions  are  continuous  across  neighboring  elements, 
and  given  that  the  problem  is  solved  on  a  closed  manifold, 
the  boundary  integral  fr  vanishes.  On  the  other  hand,  in 
the  case  of  the  DG  method,  the  boundary  integral  remains 
because  of  the  non-zero  inter-element  fluxes.  Furthermore, 
the  surface  integrals  for  DG  are  defined  element-wise  over 
Qe  rather  than  globally  over  O.  This  makes  Equations  (13) 
and  (14)  a  system  of  Ne  equations  that  are  coupled  via  the 
flux  of  Equation  (14).  Both  options  are  available  within  the 
same  code  used  for  this  study  and  both  methods  share  most 
of  the  numerical  machinery. 

The  integrals  above  are  solved  element-wise  on  the  discrete 
polyhedral  approximation  kth.  The  discrete  domain  Qh  is 
defined  by  the  union  of  Ne  non-overlapping  quadrilateral 
elements  S+.  For  every  element  we  define  a  bijective 
transformation  J-^h  :  (x,  y,  z)  — >  (£,  r],  CJ)  that  maps  the 
physical  coordinate  system,  (x,  y,  z)  to  the  reference  system 
(£,?hC)  and  is  such  that  the  reference  element  Og(f,7/)  = 
[—1,1]  x  [—1,1]  lies  on  the  spherical  surface.  £  is  thus  the 
spherical  radius  whose  discrete  values  identify  a  spherical 
shell  (e.g.  £  =  1  is  the  sphere  of  unit  radius.)  The  Jacobian 
matrix  of  the  transformation  has  components  Jl  =  d^T1 
and  determinant  |J|  =  c^-x  •  G,  where  G  =  9jx  x  dvx.  is 
the  surface  conforming  component  of  the  transformation 
(Song  and  Wolf  1999;  Giraldo  2001). 

Given  the  definition  of  the  reference  element,  the  element¬ 
wise  solution  can  be  approximated  by  the  expansion 

(N+1)2 

qN{x,t) \ne=  ^  +(+!'e1(x))'?iv(xA:,I),  e  =  l,...,lVe 

k= 1 

(15) 

where  (N  +  l)2  is  the  number  of  collocation  points  within  the 
element  of  order  N,  and  ipk  are  the  interpolation  polynomials 
evaluated  at  point  k.  The  basis  functions  ipk  are  constructed 
as  the  tensor  product  of  the  one-dimensional  functions 
/q(f(x))  and  hj(r)(yi))  as: 

ipk  =  hi(£(x))  ®  hj(r](x)), 
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Vi,  j  =  0, N.  /q(£(x))  and  hj(r](x))  are  the  basis  functions 
associated  with  the  AT  +  1  LGL  points  G  and  r)j ,  respectively, 
given  by  the  roots  of 

(1  -  £2)-Pat(£)  =  0, 

where  P'N{£)  is  the  derivative  of  the  iVtft-order  Legendre 
polynomial.  Given  these  definitions,  the  one-dimensional 
Lagrange  polynomials  /i;(£)  are 


is  diagonal  (assuming  inexact  integration),  with  an  obvious 
advantage  if  explicit  time  integration  is  used. 

Concerning  the  discontinuous  Galerkin  approximation,  the 
problem  at  hand  is  solved  only  locally,  and  unlike  the  case 
of  spectral  elements,  the  flux  integral  in  Equation  (14)  must 
be  discretized  as  well.  The  element-wise  counterpart  of  the 
matrix  problem  (18)  is  then  written  as: 

=  -(Ms-e)TF(qe)  +  (De)rF(qe)  +  S(qe),  (19) 


hi(0 


i  (i  -  mi® 

N(N+  1)  (£-6)^(6)' 


The  functions  hj(rj)  are  computed  in  the  same  way. 

The  same  expansion  above  is  used  to  construct  the 
derivatives.  We  write: 


where  we  obtain  Ms,e  =  (Me)-1Ms’e  from  the  element 
boundary  matrix,  Ms,e,  and  the  element  mass  matrix, 
Me.  Out  of  various  possible  choices  for  the  definition  of 
the  numerical  flux  F(q)  in  Equation  (19),  we  selected  the 
Rusanov  flux  that,  for  the  inviscid  part  of  the  flux,  is 


dqN{*,t) 

dt 

dqN{x.,t) 

c)x 


(N+ 1)2 


be  =  E 


dqN(xk,t) 

dt 


dlpkiPn  1(x)) 

=  E  - 5^ - qN(xk,t).  (16) 


ax 


The  integrals  are  computed  by  numerical  quadrature  on 
the  reference  element  as  follows: 


Ini 


q(x,  £)dx  =  /  q(£,t)|J(£)|d£  w 
Jfy 

N+l  N+l 

EE 


(17) 


*= i  l=i 


F(q)  =  i  [Fjv(q^)  +  Fiv(q^)  -  |A|(qg  -  q^)n]  ,  (20) 

where  |A|  =  |n  ■  u|  +  +cj)  is  the  maximum  wave  speed  of 
the  shallow  water  equations.  By  identifying  an  intrinsic 
tangential  direction  in  each  inter-element  boundary  edge,  q^ 
and  qA  denote  the  solution  in  the  left  (L)  and  right  (R) 
elements,  respectively.  For  a  detailed  description  of  the  DG 
algorithm  that  is  used  for  this  study,  the  reader  is  referred 
to  Kopera  and  Giraldo  (2014). 


where  are  the  Gaussian  quadrature  weights.  In  the 

case  of  spectral  elements,  the  substitution  of  the  expansions 
(15,16)  into  the  weak  form  (14)  yields  the  semi-discrete  (in 
space)  matrix  problem 

^=DTF(q)  +  S(q)  (18) 


4-2.  Time  integration 

The  solution  of  the  systems  of  ordinary  differential 
equations  (18,19)  is  advanced  in  time  by  a  third-order 
strong  stability  preserving  Runge-Kutta  method  (SSP- 
RK3)  (Cockburn  and  Shu  2001;  Spiteri  and  Ruuth  2002) 
that  yields  the  solution  at  time-step  n+l  as 


where,  for  the  global  mass  and  differentiation  matrices,  M 
and  D,  we  have  that  D  =  M“1D.  The  global  matrices 
are  obtained  from  their  local  counterparts,  Me  and  De, 
by  means  of  the  direct  stiffness  summation  operation  that 
maps  the  local  degrees  of  freedom  of  an  element  fig  to  the 
corresponding  global  degrees  of  freedom  in  and  adds 
the  element  values  in  the  global  system.  By  construction,  M 
©2013  Royal  Meteorological  Society 


qfc  =  agq"  +  « fq^1  +  R(qfc“1),  for  k=l,...,K, 

where  R{ q)  indicates  the  right-hand  sides  of  equations 
(18,19),  K  is  the  number  of  stages  of  the  RK  method, 
and  q°  =  q™,  qA  =  qn+1.  All  the  computations  are  executed 
with  a  maximum  Courant-Friedrichs-Lewy  number  (CFL, 
Courant  et  al.  (1967))  approximately  0.5.  This  value  ensures 
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the  stability  of  the  scheme  given  that  the  time  step  At  is 
such  that 


At  <  CFL  min  [|u  ■  x|  +  <t>s/X  •  xj_1- 
xenh 


0 


Lexp 


_ 1 _ 

(v’-yoHv-vi) 


if  P  <  Po  or  p  >  Pi 
if  ip0  <  <P  <  <Pi, 

(21) 


The  CFL  oc  Af|u|/A|x|min,  where  A|x|mirl  gets  smaller 
with  the  order  N. 

The  local  grid  distortion,  x,  is  defined  as  the  vector 
l£sl  .  \Vx\  |4|  ,  M  16*1  M' \ 

Ar  Aj/’Ar  AVAr  a?7  ) 

where  £,x,y,z  and  rix,y,z  are  the  metric  terms  and  A£,  A 77  are 
the  local  average  grid  sizes  across  the  surface. 


where  umax  =  80ms-1  is  the  maximum  zonal  velocity 
and  (pa,pi)  =  (7t/7,  7t/2  —  po)  are  the  lowest  and  highest 
latitudes  that  delimit  the  jet.  At  the  jet  mid-point,  where 
p  =  7t/4,  the  non-dimensional  parameter  en  =  exp[— 4/(</?j  — 
Po)2\  normalizes  the  jet  magnitude  to  «max.  From  the  initial 
zonal  flow  given  by  (21),  the  height  <j>  is  found  by  solving 


In  the  absence  of  dissipation  (e.g.,  artificial  diffusion),  the 
high  order  numerical  solution  of  advective  problems  may  be 
characterized  by  spurious  high  frequency  modes.  To  preserve 
full  stability,  along  with  respecting  the  CFL  condition  we 
applied  the  standard  spatial  Boyd-Vandeven  filter  (Boyd 
1996)  at  every  time  step.  The  filter  constant  is  chosen  to 
reduce  by  5%  the  highest  modes  only. 


9<t>(,<p)  =  9$  o 


/  +  u{p') 


tan(^j') 

r 


(22) 


where  the  integral  on  the  right-hand  side  is  computed  by 
numerical  quadrature.  As  in  the  reference,  po  is  chosen 
to  give  a  global  mean  layer  depth  of  10  km.  The  height 
perturbation  of  test  (ii)  is  the  bump  whose  shape  is  given 
by 


5.  Tests 


To  evaluate  the  behavior  of  the  solution  on  the  three  grid 
families,  we  first  solve  the  nonlinear  zonal  dynamics  test 
by  Galewsky  et  al.  (2004)  whose  exact  analytic  solution 
is  unknown,  and  follow  with  the  steady  state  nonlinear 
geostrophic  flow  by  Williamson  et  al.  (1992)  that,  on  the 
contrary,  admits  an  analytic  solution.  The  tests  are  run  at 
different  resolutions  using  elements  of  order  8;  based  on  the 
analysis  by  Gabersek  et  al.  (2012),  8  is  the  order  that  we  are 
likely  to  use  when  running  simulations  with  NUMA. 

The  test  by  Galewsky  et  al.  (2004)  consists  of  (i) 
an  unperturbed  mid-latitude  zonal  jet  on  top  of  a 
geostrophically  balanced  geopotential  height,  and  of  (ii)  its 
perturbed  counterpart  where  the  perturbation  is  triggered  by 
a  bump  in  the  geopotential.  In  either  case,  the  initial  zonal 
velocity  and  height  are  a  function  of  the  latitude  p  (measured 
in  radians  unless  otherwise  specified)  as 
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<t>'{ a,  p)  =  4* cos  (<p)exp[— A2/a2]exp[— (<^2  -  pi)2 / fi2} 

for  —  7r  <  A  <  7r,  (23) 

where  A  is  the  longitude,  0=12Om,  p2  =  7t/4,  a  =  1/3 
and  (3  =  1/15.  The  definition  of  this  bump  is  such  that  the 
perturbation  is  forced  to  zero  at  the  poles. 

The  flows  in  (i)  and  (ii)  are  simulated  along  a  time 
span  of  144  hours.  In  the  absence  of  any  perturbation,  the 
balanced  held  and  zonal  how  are  expected  to  maintain  their 
initial  state  indefinitely.  The  analytical  initial  held  is  used 
to  compute  the  L2  error  norms  of  the  numerical  solutions  as 
time  evolves. 

Although  its  solution  seems  trivial,  this  test  is  especially 
demanding  for  a  shallow  water  model  because  of  the  large 
sensitivity  of  the  solution  to  the  grid  resolution  and  geometry. 
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Unlike  other  tests  that  may  be  more  forgiving  in  this 
respect,  when  the  grid  is  not  sufficiently  fine,  grid  related 
oscillations  may  spoil  the  equilibrium  with  devastating 
consequences.  This  is  clearly  visible  in  Fig.  6,  where  the 
initially  unperturbed  and  balanced  zonal  flow  completely 
loses  its  initial  state  when  the  test  is  executed  on  the 
low  resolution  cubed-sphere  and  icosahedral  grids.  All  the 
runs  are  executed  at  the  equivalent  resolutions  AA  ~  A <p  = 
{0.625°,  1.25°,  2.5°,  5.0°}  and  for  a  varying  number  of  grid 
points  Np  =  (12000,  25000,  50000, 100000}.  The  reason  for 
the  two  sets  of  runs  will  be  clarified  in  Remark  1. 

Below,  the  keywords  Equilibrium  and  Perturbed  jet  are 
used  to  identify,  respectively,  problems  (i)  and  (ii). 

Remark  1:  equivalent  resolutions  A  few  words  are  in 
order  regarding  the  concept  of  equivalent  resolution  when 
grids  of  very  different  geometries  are  compared.  In  the  case 
of  spectral  elements  that  rely  on  LGL  nodes,  the  distance 
between  two  consecutive  grid  points  changes  within  the  same 
element.  Because  of  this,  the  equivalent  angular  resolution 
along  a  latitude  circle  subdivided  into  Ne>\  elements  of  order 
N  is  2k /(Net\N).  Given  this  definition,  we  add  that  two 
resolutions  are  approximately  equivalent  when,  within  the 
region  of  major  interest  (e.g.,  where  the  dynamics  is  most 
likely  to  develop  from  a  certain  initial  state),  the  size  of  the 
angular  resolution  of  two  different  grids  are  comparable. 
In  other  words,  when  we  say  that  the  resolution  is,  e.g., 
AA  =  0.625°,  it  means  that  the  mean  angular  distance 
between  two  meridians  in  the  proximity  of  the  jet  measures 
0.625°.  This  definition  makes  the  measurement  on  the 
icosahedral  grid  not  straightforward  because  of  the  irregular 
orientation  of  the  elements.  To  obviate  this  problem,  the 
error  norms  in  the  analysis  below  are  also  compared  with 
respect  to  the  total  number  of  grid  points. 

5. 1 .  Equilibrium 

Fig.  6  shows  the  vorticity  held  after  6  days  for 
the  unperturbed  (Equilibrium)  jet  problem  on  the 
three  different  conforming  grid  types  (rows)  and  Np  = 


{25000, 50000, 100000}  grid  points.  For  two  out  of  the 
three  grid  geometries,  the  effect  of  the  misalignment  is 
still  important  for  as  many  as  25000  grid  points,  which 
approximately  corresponds  to  an  equivalent  resolution  AA  = 
A  ip  =  1°  on  the  cubed-sphere.  The  solution  on  the  RLL 
grid  can  maintain  the  jet  with  approximately  80%  fewer 
grid  points  than  the  cubed-sphere  grid  and  50%  fewer  than 
the  icosahedral  grid.  This  is  not  reason  enough  to  draw 
conclusions  on  the  possible  superiority  of  one  grid  with 
respect  to  the  others;  rather,  it  should  only  suggest  that, 
unless  a  sufficiently  high  resolution  is  used  with  grids  that 
are  not  aligned  with  the  characteristic  flow  (which  is  likely 
to  be  the  case  in  realistic  simulations),  we  should  not  trust 
the  solution  with  high  fidelity. 

In  Figs.  7  and  8  we  plot  the  normalized  L2  error  norms 
of  the  geopotential  height  and  zonal  velocity  against  the 
number  of  grid  points  and,  respectively,  the  equivalent 
resolution  ( AA,  A<^)  (see  the  tabulated  values  of  all  the  errors 
plotted  in  the  paper,  in  the  Appendix).  These  errors  are 
computed  using  spectral  elements  on  conforming  grids.  Let 
us  focus  on  the  curves  with  continuous  lines  as  they  trace 
the  error  on  the  uniform  conforming  grids.  The  conclusions 
drawn  from  the  vorticity  contours  of  Fig.  6  are  confirmed 
by  these  error  estimates.  Furthermore,  from  Fig.  8  it  is 
interesting  to  see  that  the  L2  error  of  the  solution  on  the 
cubed-sphere  reaches  its  peak  value  at  Aip  =  AA  «  1.25°. 
This  leaves  very  little  room  for  resolution  options  when  the 
cubed-sphere  is  used  without  viscosity  at  low  resolution. 

The  results  reported  so  far  were  obtained  using  CG  on 
conforming  grids.  However,  we  also  ran  the  same  tests  using 
DG  on  every  grid.  As  expected,  the  errors  are  practically 
identical  to  those  obtained  by  CG;  the  solution  accuracy 
for  CG  and  DG  should  be  identical  because  the  solution  is 
relatively  smooth  (differences  should  only  emerge  between 
the  two  methods  in  the  presence  of  discontinuities).  For 
comparison,  we  report  the  values  of  the  CG  and  DG 
normalized  error  norms  in  Tables  3  and  5  in  the  Appendix. 
We  will  show  DG  results  on  non-conforming  grids  with  static 
rehnement  in  Sec.  5.3. 


©2013  Royal  Meteorological  Society 
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Figure  6.  Equilibrium ,  inviscid  CG  solution:  vorticity  field  (s— 1)  after  6  days.  The  solutions  on  the  cubed-sphere  (Hex,  top  row),  on  the  RLL 
(middle),  and  on  the  icosahedral  grids  (bottom  row)  are  plotted  from  left  to  right  using  25000,  50000,  and  100000  grid  points.  The  color  scale 
ranges  between  — 12  x  10-5  s-1  (dark  blue)  and  16  x  10-5  s-1. 


In  Fig.  9  we  plot  the  time  evolution  of  \\4>\\l2  during  6  days 
for  100000  point  grids.  After  the  gravity  wave  adjustment 
during  the  first  6  hours  (see  Galewsky  et  al.  (2004)),  the 
regularity  of  the  grid  geometry  still  plays  an  important  role  in 
the  evolution  of  the  error  even  at  a  resolution  that,  in  global 
mode,  is  typically  considered  fairly  high  in  an  operational 
setting. 

So  far,  the  viscous  term  in  the  equations  has  been  ignored 
(i.e.  5  =  0).  In  the  quest  for  diminishing  the  negative  effects 
of  the  low  resolution  on  certain  grids,  we  turned  viscosity 
on  and  recomputed  the  normalized  L2  error  norms  that 
we  report  in  Table  1.  We  compare  these  errors  against  the 
inviscid  estimates  using  12000  and  25000  grid  points  and 
point  out  the  instances  when  artificial  viscosity  is  either 
beneficial  or  harmful  to  the  error  measure.  To  interpret 
(c)  2013  Royal  Meteorological  Society 


these  values,  we  need  to  take  a  parallel  look  at  the  contour 
plots  of  Fig.  10.  In  the  case  of  RLL,  where  the  grid  is 
aligned  with  the  characteristic  flow,  it  was  shown  above 
that  the  resolution  does  not  significantly  affect  the  solution. 
The  opposite,  however,  occurs  in  the  case  of  the  other  two 
grids.  Not  surprising,  by  adding  viscosity,  what  was  a  bad 
inviscid  solution  seems  to  improve;  in  contrast,  the  solution 
that  was  well  behaved  in  the  inviscid  case  has  deteriorated. 
This  behavior  is  expected  for  two  main  reasons:  (1)  the 
viscous  and  inviscid  equations  are,  mathematically,  not  the 
same  set  of  equations,  and  physically  represent  two  different 
dynamics;  (2)  viscosity  is  so  large  that  it  is  damping  all  of  the 
modes  triggered  by  the  grid  on  the  one  hand  but,  in  the  same 
way,  damps  the  solution  also  where  it  should  not.  A  thorough 
analysis  of  the  diffusion  mechanisms  that  are  currently  used 
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Figure  7.  Equilibrium,  at  day  6,  inviscid:  Normalized  L 2  error  norm  of  (ft  (left)  and  us  (right)  against  the  number  of  grid  points.  The  error  is 
computed  with  respect  to  the  initial  condition  of  the  balanced  problem  on  every  grid.  In  the  legend,  the  symbols  1L  and  2L  indicate  that  the 
solution  is  computed  on  a  non-conforming  grid  with  1  or  2  levels  of  refinement.  In  these  two  cases,  the  solution  is  computed  using  DG,  otherwise, 
CG  is  used  instead.  See  tabulated  values  in  the  Appendix. 


Figure  8.  Equilibrium  at  day  6,  inviscid:  Normalized  L2  error  norm  of  0  (left)  and  us  (right)  against  the  resolution  measured  in  degrees.  The 
error  is  computed  with  respect  to  the  initial  condition  of  the  balanced  problem  on  every  grid.  As  in  Fig.  7,  the  symbols  1L  and  2L  in  the  legend 
indicate  that  the  solution  is  computed  using  1  or  2  levels  of  refinement  on  non-confirming  grids.  In  these  two  cases,  the  solution  is  obtained  using 
DG,  otherwise,  CG  is  used  instead.  See  tabulated  values  in  the  Appendix. 


in  atmospheric  problems  goes  beyond  the  scope  of  this 
paper.  Nevertheless,  it  should  be  emphasized  that  using 
an  arbitrary  diffusion  operator  may  not  actually  improve 
the  solution  quality  even  though  the  final  solution  may 
appear  smooth.  This  simple  result  stresses  the  importance 
of  deriving  a  proper  solution-based  stabilization  mechanism; 


we  are  working  on  this  topic  for  the  solution  of  the  fully 
compressible  Euler  equations,  and  will  report  our  results  in 
a  future  article.  For  the  time  being,  we  point  the  reader 
to  the  recent  works  by  Bao  et  al.  (2014),  Yu  et  al.  (2014), 
Marras  et  al.  (2014),  and  Guba  et  al.  (2014). 


©2013  Royal  Meteorological  Society 
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Figure  9.  Equilibrium,  inviscid:  time  evolution  of  the  normalized  L2  error  norm  of  <f)  (left)  and  us  (right)  relative  to  the  initial  condition  computed 
on  100000  grid  points.  As  for  the  figures  above,  1L  and  2L  correspond  to  the  DG  solution  on  the  non-conforming  grids.  On  the  other  hand,  the 
error  curves  RLL,  Hex  and  Ico  correspond  to  the  CG  solution  on  conforming  grids. 
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Figure  10.  Equilibrium  25000'.  vorticity  field  (s  1 )  at  day  6,  CG  only.  Top  row:  inviscid  solution.  Bottom  row:  artificial  diffusion  with 
v  —  lx  L05m2  s  1 .  Solution  on  the  RLL  (left),  cubed-sphere  (center),  and  icosahedral  (right)  grids. 


5.2.  Barotropic  instability 

Unlike  the  equilibrium  problem,  we  do  not  have  an  analytic 
solution  for  the  perturbed  jet.  Because  of  this,  the  solution 
at  the  highest  resolution  AA  =  Aip  =  0.625°  is  considered 
the  reference  to  compute  the  normalized  errors.  At  this 
resolution,  the  RLL  reference  solution  consists  of  180000 
grid  points.  Based  on  the  analysis  of  the  equilibrium  test 
and  the  high  sensitivity  to  the  resolution  in  the  case  of  the 
(c)  2013  Royal  Meteorological  Society 


cubed-sphere  and  icosahedral  grids,  we  considered  the  high 
resolution  RLL  run  to  be  the  true  solution  against  which  the 
comparison  should  be  performed.  It  is  a  choice  that  is  also 
supported  by  the  fact  that  the  difference  between  the  RLL 
and  Hex  solutions  first,  between  the  Hex  and  Ico  solutions 
then,  and  finally  between  the  solutions  on  the  RLL  and  Ico 
grids,  are  of  the  same  order  of  magnitude  (approximately 
vorticity=  O(10-5s-1);  the  plots  of  the  differences  are  not 
shown).  This  gives  us  confidence  that  the  comparison  with 
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Table  1.  Equilibrium:  viscous  vs.  inviscid  solution.  Normalized  L2 
error  norms  of  (<f>,  us,vs)  at  day  6  using  CG  with  artificial  viscosity 
(V)  v  =  1  X  105  m2  s-1  on  12000  and  25000  grid  points.  For  a  direct 
comparison,  the  inviscid  (I)  errors  computed  on  the  same  grids  are 
also  reported,  value  indicates  that  the  use  of  viscosity  induced  an 
error  increase  in  the  solution.  On  the  contrary,  value  indicates  that 
the  error  norm  decreased  by  using  viscosity. 


CG,  L2 

N  points 

12000-V 

12000-1 

25000-V 

25000-1 

RLL 

1.493e-3 

1.466e-4 

1.481e-3 

2.011e-5 

us 

1.575e-l 

2.701e-2 

1.548e-l 

5.162e-3 

Vs 

1.575e-l 

2.706e-2 

1.548e-l 

5.167e-3 

Hex 

<t> 

3.490e-3 

5.460e-3 

1.744e-3 

5.440e-3 

us 

3.171e-l 

5.059e-l 

1.814e-l 

5.458e-l 

vs 

3.172e-l 

5.059e-l 

1.815e-l 

5.458e-l 

Ico 

0 

7.002e-3 

1.263e-2 

3.284e-3 

1.057e-2 

us 

5.762e-l 

9.163e-l 

2.947e-l 

8.035e-l 

Vs 

5.762e-l 

9.163e-l 

2.947e-l 

8.035e-l 

respect  to  the  RLL  will  not  be  biased.  After  6  days,  the 
instability  has  fully  developed.  The  vorticity  computed  at 
high  resolution  is  plotted  in  Fig.  11  for  the  three  grids. 
There  is  no  visible  difference  among  the  three  solutions  with 
Aip  =  AA  0.625°.  The  solutions  are  practically  identical 
for  Np  >  25000  (see  Fig.  12).  However,  we  see  in  Fig.  13  that 
the  error  increases  on  the  cubed-sphere  and  Ico  grid  as  the 
resolution  is  decreased  below  25000  grid  points. 

Unlike  for  the  equilibrium  case,  the  error  measures  of  this 
problem  are  very  similar  for  the  three  different  grids.  This  is 
expected  because  of  the  highly  irregular  structure  of  the  jet 
which  is  no  longer  aligned  with  any  grid.  As  the  resolution 
is  coarsened,  the  error  is  sensibly  higher  than  its  equilibrium 
counterpart. 

In  the  next  Section,  we  proceed  by  testing  the 
implementation  for  non-conforming  grids  with  static  mesh 
refinement. 


5.3.  Non- conforming  grids 

To  test  the  DG  solver  on  non-conforming  grids,  we  define 
a  statically  refined  region  in  the  proximity  of  the  zonal 
wind.  The  high  resolution  static  refinement  is  positioned 
within  the  latitude  band  35°  <  ifiband  <  55°;  this  region 
covers  approximately  10°  above  and  below  the  central-line 
of  the  initial  jet.  Similar  to  the  uniform  case,  we  define  a  set 
of  different  equivalent  grid  resolutions  in  the  banded  region, 
(c)  2013  Royal  Meteorological  Society 


We  use  the  values  AA  band  =  Afband  ~  {0.625°,  1.25°,  2. 5°). 
After  fixing  the  resolution  in  this  region,  for  each  grid  we 
add  a  1-level  and  a  2-level  global  (de-)rehnement.  The  total 
number  of  grid  points  in  the  global  mesh  is  thus  reduced.  The 
grid  structure  is  such  that  each  refined  region  has  half  the 
resolution  of  its  neighbors.  If  two  refinement  levels  are  used, 
and  Ag>  =  AA  =  1.25°  characterizes  the  finest  portion  of  the 
mesh,  the  resolution  farthest  away  from  it  has  A<p  =  AA  = 
5.0°.  Based  on  the  errors  found  for  the  coarse  conforming 
grids,  we  did  not  go  beyond  this  value  anywhere  in  the 
domain.  A  total  of  6  non-conforming  grids  are  generated  and 
are  shown  in  Fig.  14. 

It  is  well  known  that  the  fundamental  advantage  of  static 
and  dynamic  grids  is  reflected  in  the  lower  computational 
cost  in  terms  of  floating  point  operations  since  the  number 
of  degrees  of  freedom  of  the  grid  is  lower  than  its  uniform 
counterpart.  The  normalized  L2  error  norms  of  ((f>,us,vs) 
are  plotted  in  Figs.  7  and  8  as  a  function  of  the  number  of 
points  given  a  certain  resolution  in  the  region  where  most 
of  the  dynamics  happen.  The  errors  are  computed  for  one 
and  two  refinement  levels.  We  observe  that  the  number  of 
grid  points  necessary  to  obtain  a  certain  error  is  reduced  by 
adopting  one  level  of  refinement  for  every  grid.  We  would 
expect  the  same  to  be  true  once  the  grid  is  further  coarsened 
using  two  refinement  levels,  however,  this  is  not  the  case 
in  this  particular  instance.  Due  to  the  non-linearity  of  the 
problem  and  the  motion  of  the  gravity  waves  across  the  whole 
domain  during  the  first  hours  of  the  simulation,  the  effect  of 
the  poorly  resolved  gravity  waves  in  the  southern  hemisphere 
is  clearly  felt  by  the  jet  resolved  at  a  hner  resolution.  Back 
to  the  6-day  error  evolution  (Fig.  9),  this  is  clearly  observed 
at  all  times.  This  issue  should  not  deter  us  from  using  more 
than  one  level  of  refinement  in  general;  it  simply  confirms  the 
well  known  fact  that  the  criterion  for  grid  adaptation  should 
be  very  well  matched  with  the  problem  at  hand  across  the 
entire  domain. 

5-4-  Dynamic  grid  refinement 

Both  static  and  dynamic  grid  rehnement  are  possible. 
Initially,  a  uniform  and  conforming  grid  is  built  before 
Prepared  using  qjrms4.cls 


Mar r as  et  al. 


15 


Figure  11.  Perturbed  jet:  vorticity  field  (s  1)  at  day  6,  CG  only.  The  flow  on  the  RLL  grid  (left)  the  cubed-sphere  (middle),  and  Ico  grid  (right) 
is  computed  with  Ac/?  «  0.625°. 
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Figure  12.  Perturbed  jet,  inviscid,  CG  only:  vorticity  field  (s  1 )  at  day  6.  Top  row:  Hex.  Middle  row:  RLL.  Bottom  row:  Ico.  Solutions  plotted  for 
25000,  50000,  and  100000  grid  points  (from  left  to  right).  The  color  scale  ranges  between  —12  X  1 0  5  s  1  (dark  blue)  and  16  X  1 0  5  s  1 . 


proceeding  to  a  first  grid  adaptation  around  the  main  jet 
identified  by  a  non-zero  vorticity.  As  the  time  evolves  and 
the  jet  breaks,  the  grid  is  adapted  according  to  a  criterion 
that,  in  the  case  of  this  work,  was  chosen  as  the  jet  condition 
such  that  the  vorticity  |A|  <  3e  —  5  s-1  (like  in  St-Cyr  et  al. 
(2008)).  The  initial  grid  configuration  and  the  dynamically 
adapted  grids  after  6  days  are  shown  in  Fig.  15.  We  can  state 


with  sufficient  confidence  that  the  adaptation  has  performed 
as  expected  for  all  the  grids.  The  jets  on  the  three  grids  are 
highly  comparable,  and  so  are  the  regions  of  refinement  and 
coarsening,  indicating  that  the  vorticity  amplitude  for  the 
three  grids  is  the  same  and  corresponds  to  |A|  <  3e  —  5  s-1. 


©2013  Royal  Meteorological  Society 
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Figure  13.  Perturbed  jet  at  day  6,  inviscid:  Normalized  L2  error  norms  of  c6  (left)  and  us  (right)  against  the  number  of  grid  points.  The  error  is 
computed  with  respect  to  the  8th  order  RLL  solution  at  0.625°,  which  corresponds  to  180000  points.  See  tabulated  values  in  the  Appendix. 


Figure  14.  High-order  (curved  elements)  non-conforming  grids  with  1-level  (top  row)  and  2-level  (bottom  row)  refinements.  RLL  (left  column), 
cubed-sphere  (middle  column),  and  Ico  (right  column).  The  highest  resolution  in  the  main  region  of  interest  is  AA  =  A  if  ~  0.625°. 


5.5.  Williamson’s  case  2 

To  further  analyze  the  grid  effects,  we  run  the  case  2  by 
Williamson  et  al.  (1992)  which  has  an  analytic  solution.  The 
test  is  run  using  a  45°  north-wise  rotation  of  the  flow 
with  respect  to  the  grid  with  Np  =  (25000,  50000, 100000} 
(c)  2013  Royal  Meteorological  Society 


points.  For  a  direct  comparison  with  St-Cyr  et  al.  (2008), 
we  also  run  this  test  using  a  uniform  resolution  of  2.5°  and, 
additionally,  of  5°.  The  error  norms  are  plotted  in  Fig.  16. 
Although  the  steady-state  solution  is  greatly  affected  by  all 
three  rotated  grids,  smaller  error  norms  are  found  in  the 
case  of  the  cubed-sphere.  Because  8  of  the  14  patches  of 
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Figure  15.  Dynamic  grid  refinement.  DG  solution.  Top  row:  initial  grid  and  vorticity  field.  Bottom  row:  adapted  grid  and  vorticity  at  day  6.  From 
left  to  right:  cubed-sphere,  RLL,  icosahedral.  The  refinement  was  triggered  by  a  vorticity  value  of  |A|  <  3e  —  5s— 1 . 


the  RLL  grid  are  built  using  a  linear  interpolation  (TFI) 
followed  by  a  projection  to  build  the  high-order  grid  LGL 
points  on  the  sphere,  a  series  of  interpolation  errors  seem  to 
cause  the  current  RLL  grid  to  lose  accuracy  in  the  interface 
patches  (Fig.  2).  To  prove  this  conjecture,  we  compute  the 
difference  between  the  analytic  and  the  numerical  solutions 
on  the  RLL  grid  and  found  that  the  larger  errors  are  indeed 
concentrated  on  the  grid  points  within  the  interface  patches 
(plot  of  the  difference  not  shown).  Since  the  flow  is  not 
aligned  in  either  one  of  the  three  grids,  we  expect  that  an 
optimized  RLL  grid  could  improve  the  solution  to,  at  least, 
the  error  norms  of  the  cubed-sphere.  Nonetheless,  the  error 
norms  are  quantitatively  comparable  to  the  results  obtained 
by  St-Cyr  et  al.  (2008). 

6.  Conclusions 

We  analyzed  the  behavior  of  a  unified  continu¬ 
ous/discontinuous  element-based  Galerkin  model  to  solve 
the  shallow  water  equations  on  six  types  of  spherical  grids. 
We  adopted  the  cubed-sphere,  a  type  of  reduced  Lat-Lon 
(RLL)  grid,  and  the  icosahedral  grid  in  conforming,  non- 
conforming,  static,  and  dynamic  configurations.  Among 
©2013  Royal  Meteorological  Society 


others,  one  advantageous  characteristic  of  element-based 
Galerkin  methods  is  that  they  are  not  constrained  by  the 
logical  structure  of  the  mesh.  For  this  reasons,  it  is  not 
surprising  that  the  spectral  element  and  discontinuous 
Galerkin  solutions  showed  almost  equivalent  error  measures 
when  executed  in  conforming  mode  on  equivalent  grids.  In 
the  case  of  non-conforming  grids,  we  only  showed  the  DG 
results  for  brevity  since  the  CG  and  DG  results  are  expected 
to  be  quite  similar  because  the  flow  analyzed  in  this  paper 
has  no  discontinuities. 

To  show  the  flexibility  of  our  element-based  Galerkin 
model  with  respect  to  the  mesh  and  see  how  the  solution 
may  be  sensitive  to  the  grid  configuration,  we  used  two 
non-linear  zonal  flow  tests.  First,  we  solved  the  time- 
dependent  geostrophic  jets  proposed  by  Galewsky  et  al. 
(2004).  Galewsky ’s  tests  consist  of  (i)  an  unperturbed  mid¬ 
latitude  zonal  jet  on  top  of  a  geostrophically  balanced 
geopotential  height,  and  (ii)  its  perturbed  counterpart  where 
the  perturbation  is  triggered  by  a  bump  in  the  geopotential. 
As  suggested  by  Galewsky  et  al.  (2004),  both  the  inviscicl 
and  viscous  solutions  with  artificial  diffusion  were  computed. 
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Figure  16.  Williamson’s  case  2.  Normalized  L2  and  Loo  errors  of  c6  at  day  5  using  CG.  Left:  error  against  the  total  number  of  grid  points.  Right: 
error  against  the  resolution.  See  tabulated  values  in  the  Appendix. 


Because  no  exact  solution  exists  for  Galewky’s  tests,  we 
also  tested  the  code  against  the  steady-state  non-linear 
geostrophic  flow  of  Williamson  et  al.  (1992)  which  admits 
an  analytic  solution. 

Although  we  have  shown  that  our  unified  element-based 
Galerkin  shallow  water  model  can  handle  any  computational 
grid,  which  implies  that  highly  optimized  arbitrary  grids 
could  be  generated  offline  and  passed  onto  the  code  without 
additional  effort,  we  also  found  that  the  large  problem- 
dependence  of  the  solutions  suggests  that  the  users  of  any 
atmospheric  model  should  be  able  to  select  the  grid  based 
on  their  knowledge  of  the  flow  to  be  solved. 

Possibly,  dynamic  grid  adaptivity  is  the  real  solution  to 
always  obtain  a  locally  optimal  grid  configuration. 
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A.  Error  tables 

In  this  Appendix,  we  report  the  tabulated  values  of  the  error 
curves  plotted  in  Figs.  7,  8,  13  for  Galewsky’s  tests,  and  Fig. 
16  for  Williamson’s  test  case  2. 
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Table  4.  Equilibrium.  As  table  3,  but  on  the  Ico  grid.  Because  of  the  geometrical  structure  of  Ico,  the  correspondence  between  the  equivalent 
resolution  and  the  number  of  subdivision  indicated  by  N6 , . . . ,  N1  is  not  as  straightforward  as  it  is  for  Hex  and  RLL. 


Geometry 

Approximate  AA 

N6 

0.625° 

N5 

N4 

1.25° 

NS 

N2 

2.5° 

N1 

5.0° 

0 

5.064e-6 

2.299e-5 

2.593e-4 

6.157e-3 

1.261e-2 

1.033e-2 

Us 

4.381e-4 

2.039e-3 

2.313e-2 

5.169e-l 

9.166e-l 

7.774e-l 

Vs 

4.381e-4 

2.039e-3 

2.313e-2 

5.169e-l 

9.166e-l 

7.774e-l 

Table  2.  Equilibrium.  Normalized  L2  error  of  (< f>,us,vs )  at  day  6 
using  CG.  The  error  is  computed  with  respect  to  the  analytical 
initial  fields  given  by  equations  (21)  and  (22),  and  is  reported  against 
the  total  number  of  grid  points. 


Table  6.  Equilibrium  on  non-conforming  grids  with  1-level 
refinement  (DG  only):  normalized  L2  error  of  (<j), us,vs)  at  day  6. 
The  error  is  computed  with  respect  to  the  analytical  initial  fields 
given  by  Eqs.  (21,22). 


CG,  L2 

N  points 

100000 

50000 

25000 

12000 

RLL 

0 

4.878e-8 

1.802e-7 

2.011e-5 

1.466e-4 

us 

5.910e-5 

1.586e-4 

5.162e-3 

2.701e-2 

Vs 

5.910e-5 

1.586e-4 

5.167e-3 

2.706e-2 

Hex 

0 

6.099e-5 

4.896e-4 

5.440e-3 

5.460e-3 

Us 

6.679e-3 

5.346e-2 

5.458e-l 

5.059e-l 

Vs 

6.679e-3 

5.346e-2 

5.458e-l 

5.059e-l 

Ico 

0 

2.230e-5 

2.593e-4 

1.057e-2 

1.263e-2 

Us 

2.040e-3 

2.313e-2 

8.035e-l 

9.163e-l 

Vs 

2.040e-3 

2.313e-2 

8.035e-l 

9.163e-l 

Table  3.  Equilibrium.  As  table  2,  but  the  error  is  reported  against 
the  equivalent  resolution  A  ip  =  A  A. 


DG  1-level,  L2 

RLL-NC 

0.625/1.25° 

1.25/2.5° 

2. 5/5.0° 

N  points 

43000 

13000 

5000 

0 

2.207e-7 

1.548e-5 

1.337e-3 

Us 

1.647e-4 

3.960e-3 

1.734e-l 

Vs 

1.647e-4 

3.960e-3 

1.730e-l 

Hex-NC 

0.625/1.25° 

1.25/2.5° 

2. 5/5.0° 

N  points 

31000 

10000 

1500 

0 

1.680e-4 

5.210e-3 

8.513e-3 

Us 

1.844e-2 

5.372e-l 

5.783e-l 

Vs 

1.844e-2 

5.372e-l 

5.783e-l 

Ico-NC 

0.625/1.25° 

1.25/2.5° 

2. 5/5.0° 

N  points 

61000 

15000 

4000 

0 

8.881e-7 

3.241e-3 

1.348e-2 

Us 

8.294e-5 

2.879e-l 

9.729e-l 

Vs 

8.294e-5 

2.879e-l 

9.729e-l 

Table  7.  As  Table  6,  but  with  2-level  refinement  (DG  only). 


CG,  L2 


RLL 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

180000 

49000 

16000 

3000 

0 

1.377e-8 

5.017e-7 

1.106e-4 

2.056e-3 

Us 

2.185e-5 

3.167e-4 

2.553e-2 

4.694e-l 

Vs 

2.185e-5 

3.167e-4 

2.555e-2 

4.694e-l 

Hex 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

124000 

31000 

9600 

1500 

0 

5.109e-5 

4.721e-3 

6.336e-3 

1.038e-2 

us 

5.601e-3 

4.709e-l 

4.832e-l 

6.918e-l 

vs 

5.601e-3 

4.709e-l 

4.832e-l 

6.918e-l 

Table  5. 

Equilibrium.  As  in  table  3,  but  using 

DG. 

DG,  L2 


RLL 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

180000 

49000 

16000 

3000 

0 

1.384e-8 

6.556e-7 

1.203e-4 

2.151e-3 

us 

2.189e-5 

3.261e-4 

2.636e-2 

4.794e-l 

Vs 

2.189e-5 

3.260e-4 

2.640e-2 

4.794e-l 

Hex 

0.625° 

1.25° 

2.5° 

5.0° 

N  points 

124000 

31000 

9600 

1500 

0 

6.649e-5 

4.948e-3 

6.397e-3 

1.540e-2 

us 

7.291e-3 

4.933e-l 

4.972e-l 

9.369e-l 

Vs 

7.291e-3 

4.933e-l 

4.972e-l 

9.369e-l 

Conf.  Spectral  and  High  Order  Methods.  Houston  Journal  of 
Mathematics,  Houston,  Texas,  pp.  267-276. 


DG  2-level,  L2 

RLL-NC 

0.625/1.25/2.5° 

1.25/2.5/5.0° 

N  points 

13000 

5000 

<t> 

1.265e-5 

9.183e-4 

Us 

3.964e-3 

1.599e-l 

Vs 

3.964e-3 

1.599e-l 

Hex-NC 

0.625/1.25/2.5° 

1.25/2.5/5.0° 

N  points 

9600 

1500 

0 

5.187e-3 

8.398e-3 

Us 

5.364e-l 

6.348e-l 

Vs 

5.364e-l 

6.348e-l 

Ico-NC 

0.625/1.25/2.5° 

1.25/2.5/5.0° 

N  points 

15000 

4000 

0 

3.670e-3 

1.298e-2 

Us 

3.249e-l 

9.695e-l 

Vs 

3.249e-l 

9.695e-l 

Table  8.  Williamson’s  case  2.  Normalized  L2  and  Loo  errors  of  (j>  at 
day  5  using  CG. 


N  points 

100000 

50000 

25000 

L2 

Hex  (j> 

6.821e-13 

1.434e-12 

2.650e-ll 

RLL  <j> 

3.631e-8 

3.744e-8 

1.897e-7 

Ico  4> 

4.997e-8 

1.051e-7 

1.415e-7 

Lqq 

Hex  (j) 

3.634e-13 

2.359e-12 

5.575e-ll 

RLL  (j) 

7.693e-10 

1.398e-9 

2.721e-8 

Ico  4> 

4.857e-9 

4.633e-9 

2.238e-8 

Browning  GL,  Hack  JJ,  Swarztrauber  PN.  1989.  A  comparison  of 


three  numerical  methods  for  solving  differential  equations  on  the  Cockburn  B,  Shu  CW.  2001.  Runge-kutta  discontinuous  Galerkin 
sphere.  Mon.  Wea.  Rev.  117:  1058-1075.  methods  for  convection-dominated  problems.  J.  of  Sci. 
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Table  9.  Williamson’s  case  2.  Normalized  L2  and  Loo  errors  of  <j>  at 
day  5  using  CG  with  uniform  resolutions  of  2.5°  and  5.0°. 


^2 

-^OO 

2.5° 

Hex  (j> 

8.109e-10 

1.6353e-9 

RLL  0 

2.081e-7 

1.056e-9 

5° 

Hex  <j> 

2.508e-7 

1.348e-7 

RLL  0 

1.102e-5 

2.311e-6 
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